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Abstract. The coupled Dyson-Schwinger equations for the gluon and ghost 
propagators are investigated in the Landau gauge using a two-loop improved 
truncation that preserves the multiplicative renormalizability of the propaga- 
tors. In this truncation all diagrams contribute to the leading order infrared 
analysis. The infrared contributions of the nonperturbative two-loop diagrams 
to the gluon vacuum polarization are computed analytically, and this reveals 
that infrared power behaved propagator solutions only exist when the squint 
diagram contribution is taken into account. For small momenta the gluon and 
ghost dressing functions behave respectively like (p^)^^ and (p^)"'^, and the 
running coupling exhibits a fixed point. The values of the infrared exponent 
and fixed point depend on the precise details of the truncation. The coupled 
ghost-gluon system is solved numerically for all momenta, and the solutions 
have infrared behaviors consistent with the predictions of the infrared anal- 
ysis. For truncation parameters chosen such that k = 0.5, the two-loop im- 
proved truncation is able to produce solutions for the propagators and running 
coupling which are in very good agreement with recent lattice simulations. 



1 Introduction 

In the Standard Model of the strong, weak and electromagnetic forces, the inter- 
actions are quantitatively described by gauge field theories. Quantum chromody- 
namics is a non-Abelian SU(3) gauge theory, and the proof of its renormalizability 
[1] and discovery of ultraviolet asymptotic freedom [2] have been milestones in its 
acceptance as theory of the strong interaction. For large momenta, perturbation 
theory seems an appropriate calculational tool, as the coupling becomes small. 
However, for small momenta the coupling grows large and adequate methods 
have to be used to study nonperturbative phenomena like confinement, chiral 
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symmetry breaking and fermion mass generation [3]. One such method is the 
study of the Dyson- Schwinger equations (DSE) [4], and their phenomenologi- 
cal apphcations to hadronic physics is a subject of growing interest [5, 6]. The 
gluon self-interaction and existence of ghost fields are remarkable features of 
non-Abelian gauge theories, and the study of the infrared behavior of the gluon 
and ghost propagators, and of the running coupling has become a major topic 
in the quest to understand the nonperturbative phenomena in QCD. 

Early studies of the Dyson-Schwinger equation for the gluon propagator in 
the Landau gauge concluded that the gluon propagator is highly singular in the 
infrared [7, 8, 9, 10]. However, these studies neglected any contribution of the 
ghost fields, and posited cancellations of certain leading terms in the equations. 
It is therefore far from certain that these solutions possess the correct infrared 
behavior of QCD. Nevertheless, these solutions are at the origin of successful 
phenomenological applications, which can be understood by observing that they 
generate the necessary integration strength in the kernels of the gap and Bethe- 
Salpetcr equations [11]. 

The important role played by the coupling of ghost and gluon fields to gener- 
ate a consistent infrared behavior of QCD in the Landau gauge was first demon- 
strated in a study of the coupled propagator equations [12]. The propagator 
dressing functions were found to behave like powers for small momenta: the 
gluon as (p^)^" and the ghost as {p'^)~'^. Furthermore, the strong running cou- 
pling is related to both propagators and exhibits an infrared fixed point ao whose 
value depends on k. Initial studies [12, 13] truncated the system using Ansdtze 
for the ghost-gluon and threc-gluon vertices, neglecting the two-loop gluon vac- 
uum polarization diagrams and introducing an angular averaging procedure. The 
infrared analysis was performed analytically and yielded k = 0.92, ao = 9.5 [12], 
and K = 0.77, ao = 11.5 [13], depending on the details of the truncation. The cou- 
pled equations were also solved numerically for all momenta, and the numerical 
solutions have the infrared behavior predicted by the analytic infrared analysis. 
The uniqueness of the solutions of ref. [13] was recently confirmed through an 
extended analysis of the infrared and ultraviolet asymptotic behaviors of the 
truncation [14, 15]. A discussion of the Faddeev-Popov quantization at the non- 
perturbative level and of the Gribov horizon problem can be found in ref. [16]. 
Comparison with lattice simulations [17, 18, 19] is quite encouraging as these 
studies have revealed quite similar infrared behaviors for the propagators, i.e. a 
vanishing gluon dressing function and a singular ghost dressing function. Never- 
theless, quantitative discrepancies in the value of k, fixed point of the coupling 
and behavior of the gluon propagator at intermediate momentum hint to the 
need of improved truncations in the Dyson-Schwinger equation studies. 

As was shown in ref. [20] the recipe used for angular averaging influences the 
value of the exponent n and even the existence of infrared power law solutions, 
and a method was developed therein to treat the angular integrals exactly. The 
infrared analysis revealed an exponent k = 1, agreeing with the Gribov form of 
the gluon and ghost propagators [21], and a smaller fixed point value ao = 4.2. 
Unfortunately, the numerical analysis could not find solutions satisfying the in- 
tegral equations over the whole momentum range, and the question whether 
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this infrared behavior is indeed a property of QCD or a spurious solution of the 
infrared analysis was left unanswered. Clearly, the existence of power laws satis- 
fying the infrared consistency conditions is not sufficient to ensure the existence 
of an overall solution of the integral equations. 

To avoid spurious unrenormalizable ultraviolet quadratic divergences, refs. 
[12, 13, 20] extracted the gluon dressing function by contracting the gluon equa- 
tion with the Brown-Pennington (BP) tensor [10]. In more recent work [22, 23] 
the contraction was performed using an arbitrary tensor, and the infrared anal- 
ysis yielded an additional solution for n slightly larger than 0.5 when using a 
non-BP tensor. These studies were restricted to an infrared analysis, and no 
full numerical solution was given. In fact, a straightforward analysis shows that 
these truncations are plagued by unrenormalizable ultraviolet quadratic diver- 
gences, and, as will be shown later in this paper, the new infrared solution for k 
is spurious as its existence is directly related to that of the ultraviolet quadratic 
divergences. 

This truncation was later enhanced [24] by removing the ultraviolet quadratic 
divergences using a prescription specially constructed to leave the infrared analy- 
sis of ref. [22] unchanged, thus retaining the spurious k, ^ 0.5 solution. The merit 
of this study is to produce the first numerical solution to the integral equations 
with full angular integrals. However, as will be discussed later, its prescription 
to remove the quadratic divergences is not satisfying, and subtracting these in a 
consistent way removes the spurious infrared solution and destroys the existence 
of the numerical solution. 

For all the above-mentioned truncations, following in the footsteps of ref. [12], 
only the ghost loop contributes to the leading infrared power of the gluon DSE. At 
the same time, all these truncations violate the multiplicative renormalizability of 
the propagators. In ref. [25] I proposed a novel truncation of the Dyson- Schwinger 
equations that respects the multiplicative renormalizability of the propagators 
and satisfies their resummed leading order perturbative behavior. A major new 
feature of the truncation is that all vacuum polarization diagrams of the gluon 
DSE contribute to its leading infrared power. An analytic infrared analysis of 
the coupled ghost-gluon system was performed taking into account the ghost and 
gluon loops in the gluon vacuum polarization, but still neglecting the two-loop 
diagrams. The study could not find infrared power laws simultaneously solving 
the gluon and ghost equations. Although the two-loop diagrams also contribute 
to the leading infrared power, their numerical contribution could not yet be 
evaluated at that time, and it was suggested that a more detailed study of the 
missing two-loop contributions is needed to get a better understanding of the 
infrared behavior of QCD. 

In the present paper the MR truncation of ref. [25] for the coupled ghost- 
gluon DSEs in pure Yang-Mills theory is generalized further. The various infrared 
contributions to the ghost self-energy and gluon vacuum polarization are com- 
puted analytically. This includes an analytic solution for the two-loop infrared 
contributions of the sunset and squint diagrams. The contribution of the squint 
diagram proves to be such that propagator power laws satisfy the infrared con- 
sistency conditions imposed by ghost and gluon equations, and the values of the 
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infrared exponent and fixed point of the coupling depend on the details of the 
truncation. 

For the numerical study, the two-loop integrals are replaced by effective one- 
loop integrals with same infrared behavior and consistent subleading ultraviolet 
behavior. The equations are then solved numerically for all momenta, and the 
numerical solutions agree with the predicted infrared and ultraviolet behaviors. 
The results for ghost and gluon propagators, and for the running coupling are 
also compared to recent results from lattice simulations [19]. 



2 The Equations 

The Dyson- Schwinger equations for the ghost and gluon propagators can be 
diagrammatically written as 



r 




q 




with the labelled momenta flowing from left to right. Several contributions to the 
gluon vacuum polarization were neglected in previous studies, and here the scope 
of the investigation is extended by considering its complete form including the 
contributions from the ghost loop, gluon loop, three- gluon loop and four-gluon 
loop^. Written out explicitly, the Euclidean space equations are 



A«''(p) 



-1 




-1 







-1 



-n"t(p), 



(3) 

(4) 

(5) 



with ghost self-energy integral^ 


^The tadpole contribution is independent of the external momentum, and vanishes after sub- 
traction of the renormalized equations. 
7 = /A/(27rr 
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and gluon vacuum polarization 

K(P) = n'r' (p) + n^^/^p) + u%'-\p) + uffj-\p) , (6) 

where 

n^r'W = (-1) / G-^"H-r,g)A«i^i(9)G°'''^'^(g,-0A''^"^(-r), (7) 

Jg 

^lr\p) = \ j I rXZSr^iP, -13, -12, -Qi)DZXin)DZ%i<l2) 

" -Jqi Jq2 (9) 

X Ktk^l2)D^^X^q,)TlfJX^^^ . (10) 

Here D is the gluon propagator, A the ghost propagator, T^f the three-gluon 
vertex, F^^ the four-gluon vertex, G the ghost-gluon vertex, the zero superscript 
denotes bare quantities, and the definitions of the momenta can be read off from 
the diagrammatic Eqs. (1) and (2). The extra factor (—1) in front of the ghost 
loop (7) is due to the Fermi statistics of the ghost field. The momenta in the 
three- and four-gluon vertices arc incoming momenta, the ghost momentum in 
the ghost-gluon vertex is incoming, and the anti-ghost momentum is outgoing. 
The integrals are regulated using an ultraviolet cutoff A. 

The general expressions for the full gluon and ghost propagators in a covariant 
gauge ^ can be written as 



^2 1 ^2 """^ „4 



p. J p. p-. 



(11) 



Dt{p) 

A"^(p) = (12) 

pA 



^2\ 



where the scalar functions F{p^) and G{p'^) are the gluon and ghost dressing 
functions, and the bare propagators correspond to F{p^) = 1 and G{p^) = 1. 
The full vertices can be decomposed in their color structures: 

Gf'{p,p') = igor'%G^^{p,p') , (13) 
^%f%p, q, r) = -igor'^rXiP, g, r) , (14) 

r^^.';r(P, Q, r, s) = -giir'^f^r^^AP, Q, r, s) + f"^f'X%{p, r, q, s) 

where go is the bare coupling and p' is the anti-ghost momentum flowing out of 
the ghost-gluon vertex. After separation of the color structure, the tensor F^^ is 
antisymmetric under exchange of any two arguments, while F^^ is antisymmetric 
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under exchange of 1 <-> 2 or 3 <-> 4, and symmetric under exchange of (1,2) 
(3,4). The bare vertices correspond to 

G%{p,p') = 6x^, (16) 
^fiup'^P^ 9' ^) = (P- 9)p V + (9 - r)i_t6i,p + (r - p)Jpij, , (17) 

^t^'paiP^ 9, r, s) = 6p,p6ya - ^ixa^up , (18) 

in Eqs. (13), (14) and (15). 

The Greens functions are potentiahy divergent as the ultraviolet cutoff is 
taken to infinity, and therefore renormahzed fields and couplings are introduced. 
Proceeding as described in ref. [25], the exact equations for the renormalized 
dressing functions in the Landau gauge can be written as 

^ Z3(M^A2)-S(p2), (19) 



^ Z3(/x2,A2)-n(p2), (20) 



where the renormalized and unrenormalized dressing functions are related 
through multiplicative renormalizability: 



'^^'^"W;^) ' ^3(M,A)-^^(^2,^2, ' (21) 



and Z3 and are the gluon and ghost field rcnormalization constants fixed by 
the renormalization conditions Fji{iJ? , ij?) = Gr{ij?,ii'^) = 1. The scalar gluon 
equation (20) was extracted from the tensor equation (4) by contracting with 
the tensor 

Vf,u{p) = 5^.-n^-^. (22) 

Throughout this paper the Brown-Pennington tensor (n = 4) will be used, un- 
less stated otherwise, such that spurious, unrenormalizable quadratic ultraviolet 
divergences are avoided [10]. The ghost self-energy contribution in Eq. (19) is 
given by 

Hp") = ^ / T{p\q\r^)G{q^)F{r^), (23) 

Z^Z^ Jq 

and the gluon vacuum polarization in Eq. (20) by 



n(p2) = n2^(p2) + n25(/) + n35(/) + u^'^ip'') , (24) 



with 



^(/) = ^ / M(p^ ,^ r^) G{q')G{r') , (25) 

Zg Jq 

2 



tf^(/) = ^ f Q{p\q\r')F{q^)F{r'), 

Z3 Jq 



(26) 



^3 



n'^(/) = ^/ / Q3{p,<li,<l2,q3)F{ql)F{ql)F{ql), (27) 

'^3^3 Jqi Jq2 
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U^'ip') = ^ / / Q4(P, 91, 92, 93, 94) F(g?)F(9i)F(gi)F(g|) , (28) 

^^3^3 Jqi Jq2 

where the contractions over color indices were performed using the formulae of 
Appendix A, and is the number of colors. The renormalized coupling was 
introduced through the renormalization of the ghost-gluon interaction 

9n = 9if^) = — , 2 A2^ 3oi^ ) ' 

where Zi is the ghost-gluon vertex renormalization constant. Equations (23)-(28) 
make use of the universality of the strong coupling stating that the renormalized 
coupling is identical for the ghost-gluon, three-gluon and four-gluon interaction. 
Furthermore, Zi = 1 in the Landau gauge, as shown by Taylor [26]. In the 
Landau gauge the various kernels in Eqs. (23)- (28) are defined by 

T{p\q\r'') = -^^gj^^{r)G^{p,q)Gl{q,p), (30) 

Mlp^q^r-") = ^ V^,,{p)G^{-r,q)Gl{q,-r), (31) 

Q{p\q\ r^) = --^^^'PAp) 5ai/3i il) 9a,/3, (r) ^^^^ 
X ^Xa, {P, -r,-q) r^|;% {-P,q,r), 

Q3(P,91,92,93) = . \ 2^2 ^M^(y) (9l)g4/32 (g2)g«3fe (gs) 

oop 9i9293 

X [(2r^aia2a3(P' "^l, "92, -93)r5'°/32/33 ("P, 9l, 92, 93) 

- ^Xia2C:iiP, -ft, -92, -q3)K'fS2P3fSi ^'P^ ^2, 93, 9l) ^^^^ 

- r^aiaaaab' "ft' "^2, -93)r^^;%^2 (-p, 93, 9l, 92)) 
-I- cyclic perm. (1, 2, 3)] , 

Q4(P,91,92,93,94) = o I 2 0^2 Pl^Ap)9ai/3i (g^gazfe (^2)34/33 (gs) 



12r9f92939i 

X 9p^aA-l4^)^%4aiiP^ -94, -9l)rJ^a3a2 ' "^S, "92) 
^ [02ft/34(-^' 92, 93, 94) + r5'3%^2(-p, 93, 94, 92) 
-2rSfftft(-P,94,92,93)], 



(34) 



where 5^1/(9) = ff^jy — q^qi^/q^- The kernels (30) to (34) contain full vertices 
which depend on the propagator functions F and G by the intermediary of the 
three-point Dyson-Schwinger equations. 

As can be seen from Eqs. (23) to (28), the renormalization point /j, only enters 
the right hand sides of Eqs. (19) and (20) through the renormalized coupling 
and the renormalization constants Z3 and Z3 . All vertices and dressing functions 
inside the integration kernels were left unrenormalized on purpose. The reason 
for this will become clear when the MR truncation is constructed in the next 
section. 
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The knowledge of the propagator dressing functions is phenomenologicahy 
important because of its relation to the renormalization group invariant running 
coupling. Prom the definitions (21), (29) of the renormalized quantities it is 
straightforward to show that in the Landau gauge the product 

a(g2,/x2) = «(/x2)Fi^(g^Ai')G|(g^M') (35) 

is renormalization group invariant, i.e. independent of lu!^ [7, 12, 13]. Hence, 
the product a{q'^,iJ?) is nothing else but the nonperturbative running coupling 
0(0^), as is easily demonstrated by choosing fi'^ = and noting that Fji{q'^, q^) = 
GR{q^,q^) = 1. The running coupling and its full nonperturbative behavior will 
be a major focus in the remaining of this paper. 



3 The Truncation 



As proposed in ref. [25] the renormalization constants and in Eqs. (23) to 
(28) are eliminated using their definitions (21), where the ratios of unrenormal- 
ized to renormalized dressing functions are independent of the chosen momentum 
p for a multiplicatively renormalizable theory. This arbitrariness in the choice of 
p will be used to construct a truncation which not only satisfies the principles 
of multiplicative renormalizability, but also satisfies the resummed leading order 
perturbative results. 

After elimination of Z3 and Z3 using Eq. (21) the full quantum corrections 
(23), (25)-(28) can be rewritten as 



n2^(p2) 
tf«(/) 



NcqI 



G{q^)F{r^) 
GMaE]FMbs] 
G(g^)G(r^) 

F{q^)F{r'^) 



G2[a2,]F0[b2,] 

F{ql)F{ql)F{ql 



T{p\q'y) 
-M{p\q\r^) 

Q{p\q\r^) 



G^[aE]i^A[bs], 
G|^[a2c]F^[b2c] 
G|[a2,]F0[b2,] 



(36) 
(37) 
(38) 



n^^(p2) = iv24 



91 Jqi 



G4[a3,]Fi[b3,] 

F(g^)F(gi)F(gi)F(9| 
G4[a4,]Fi[b4,] 



Qzkp-, qi,q2,q3 



G|j[a3,]FA[b3,] , 

(39) 

Q4{p, qi,q2,q3,q4) 



X G%MFk[Ug] , 



(40) 



where the symbolic notation 

n 



i=l 



m , 



(41) 



was introduced, which represents a product of powers of dressing functions eval- 
uated at the internal momenta. Here X stands for any of the dressing functions 
(renormalized or unrenormalized), n is the number of internal legs in the dia- 
gram, qi enumerates the momenta of the internal legs, v = (vi . . . -y^) is a vector 
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of exponents, and the superscript m gives the total sum of the exponents. For 
example, the notation G^[a3g] in Eq. (39) can be written out explicitly as 

G4[a3,] = G«3..i(g2) aa3,,2(^q2^ G«3..3(g|), 

where the exponents add up as a^g,! + 039,2 + CL3g,3 = 4. Note that in the single 
loops (36) to (38) the internal momenta qi and q2 were respectively labelled as 
q and r. 

Equations (36) to (40) are still exact, and have the peculiarity that the 
renormalization point only enters through the renormalized coupling and the 
renormalized dressing functions, as all the renormalization constants have been 
eliminated. Furthermore, the objects between square brackets contain regular- 
ized, but unrenormalized, quantities which depend on the ultraviolet regulator 
A, the external momentum and the loop momenta, but are independent of the 
renormalization point /x. 

As was discussed in detail in ref. [25] the shapes of Eqs. (36) to (40) are such 
that a truncation respecting the multiplicative rcnormalizability of the propa- 
gators and satisfying the leading order resummed perturbative results can be 
constructed in a straightforward way. The truncation consists in replacing each 
expression between square brackets in Eqs. (36) to (40) by its tree level analog 

E{p') = N^gl [ n{p',q',r')G],[e.^]Fk[h^] , (42) 

Jq 

ti'%p'') = N,gl [ Mo{p^q^r^)Gl[^i2c\F^[h2c\, (43) 

Jq 

U''^p^) = N,gl [ Qo{p'',q^r^)GU^2g]F^[h2g], (44) 

Jq 

n^^p') = N^gf, [ [ Q'i{p,qi,q2,q3)Gi[^sg]FUh3g], (45) 

Jqi Jq2 

Il^3^p')=Nlgl f [ Ql{p,qi,q2,q3,q4)G%[^^g]Fj,[h,g]. (46) 

J qi Jq2 

Because of the introduction of the product form (41) the truncation is a some- 
what generalized version of that proposed in ref. [25], where specific sets of values 
for the a's and b's of Eqs. (42)-(46) were chosen. Note that choice of the sets 
a and b can be done independently in the different self-energy and vacuum po- 
larization diagrams, as is suggested by their different subscripts. The proof that 
the truncation respects multiplicative renormalizability and the leading order 
resummed perturbative results can be found in ref. [25]. To ensure that the lead- 
ing order perturbative results are reproduced it is crucial to choose internal loop 
momenta when using Eq. (21) to eliminate the renormalization constants. This 
was provided for in the definition (41) of X™[v]. Choosing a fixed momentum, 
e.g. the renormalization point /x, would fail to reproduce the perturbative results. 

The freedom of choice of the sets a and b means that even this simple trun- 
cation procedure is far from unique. For most diagrams, concerns of naturalness 
and symmetry will fix the choice in the further treatment. For other diagrams, 
the sets will be used as parameters in the study of the solutions of the coupled 
DSEs. 
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Note that the approach developed herein is quite different from that of refs. 
[12, 23] where various expUcit corrections to the ghost-gluon and triple-gluon ver- 
tices were studied. Herein I consider effective kernel corrections to the equations, 
which guarantee that the kernels satisfy properties of MR and PT. 

The tree-level kernels of the one-loop diagrams in Eqs. (42)- (44) are 

To{x,y,z) 
Mo{x,y,z) 
Qo{x,y,z) 



where x = p"^, y = and z = r^. The two-loop kernels Q3 and Ql of Eqs. (45) 
and (46) are calculated in a straightforward way, but are more elaborate and 
therefore not explicitly given herein. 

The truncation is mainly motivated by its simplicity to achieve the two im- 
portant properties that are multiplicative renormalizability and perturbative be- 
havior. Nevertheless one would like to find a more theoretical justification for the 
truncation, which assumes that integrals with kernels containing specific combi- 
nations of full and bare vertices multiplied by various dressing functions can be 
well approximated by the same integrals over the tree-level kernels. 

Prom the Slavnov-Taylor identities (STI) for the three-gluon, four-gluon and 
ghost-gluon vertices one observes that the vertices receive corrections which are 
related to the gluon and ghost dressing functions [27, 28, 12]. That this is not only 
true for the longitudinal parts of the vertices, constrained by these identities, but 
for the complete vertex can easily be seen from the principles of multiplicative 
renormalizability. The three-gluon vertex receives corrections proportional to 
G/F, the four-gluon vertex to G'^/F and the ghost-gluon vertex to G/G, and 
the divergences of the vertices are given by the divergences of these ratios. 

All propagators and vertices in the Dyson-Schwinger equations (1) and (2) 
are full, except for one vertex per diagram which is bare to avoid double counting 
in the integral equation. When briefly assuming factorizable Ansdtze^ for the full 
vertices, a simple counting of dressing function corrections in each DSE diagram 
reveals that the MR truncation introduces an additional correction equivalent to 
a dressing of the remaining bare vertex. Hence, the MR truncation could wrongly 
be interpreted as dressing all the vertices in the nonperturbative diagrams. That 
this is not true is easily understood from the following ab absurdo argumentation. 
Placing full vertices on all interaction points would introduce double counting, 
and the perturbative results would not be reproduced. However, the MR trun- 
cation does reproduce the leading order perturbative result, hence refuting the 
allegation. As a matter of fact, the additional correction exactly ensures that 

^The factorizable Ansatz replaces the full vertex by a product of bare vertex and dressing 
functions as suggested by the STIs, which correctly reproduces its divergence. 
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the divergence of the integral matches that of the renormahzation constant in 
the seed of the integral equations, thus yielding a finite renormalized dressing 
functions as required by the left hand sides of Eqs. (19) and (20). 

It is known that the use of factorizable vertex Ansatze in the DSE does 
not reproduce the correct divergence and multiplicative renormalizability of the 
propagators. The exact dressed vertex has a structure that is much more com- 
plicated than just a multiplication of the bare vertex with ratios of dressing 
functions. Usually several tensor structures are involved, each having its own 
correction given by similar ratios of propagator dressings. This is well illustrated 
by the fermion-photon Curtis-Pennington (CP) vertex in quenched QED [29], 
constructed to satisfy the Ward-Takahashi identity and multiplicative renormal- 
izability. The novel MR truncation assumes that the integrals satisfy specific 
cancellations after contractions and integrations. That this is not as far fetched 
as it at first seems can be seen from quenched QED. Substituting the quite 
elaborate CP-vertex in the quark DSE yields kernels of great simplicity [30, 31] 
which arc of the form suggested above. Furthermore, in Appendix B an analysis 
of the DSE for the massless quark propagator in quenched QED shows how the 
one-loop corrections to the various vertex structures contribute to the self-energy 
integral and combine to produce an MR dressing function, and how the same 
result can be obtained with a modified kernel which only uses the bare vertex 
structure. 

Although the validity of the Ansatz could not yet be proven for QCD, the 
simplicity of the truncated equations, which satisfy multiplicative renormaliz- 
ability, rcsummed leading order perturbation theory and take into account the 
two-loop nonperturbative diagrams, is a good enough reason to further study 
their properties. 

Now that the truncation scheme has been defined, the search for solutions of 
the coupled integral equations can be initiated. First, their infrared asymptotic 
behavior will be investigated analytically. 

4 Infrared Analysis 

As shown in ref. [25] a major feature of the MR truncation is that all diagrams 
in the gluon vacuum polarization in Eq. (2) contribute to the leading infrared 
behavior if the propagators obey infrared power laws. Indeed, in the infrared the 
propagator power laws combine to a fixed point for the running coupling, and 
as all diagrams in the MR truncation have kernels with similar dressing function 
structure up to an additional factor of a(g^) they will all be infrared leading. This 
contrasts with the well-known fact that the sunset and squint diagrams are sub- 
leading for the ultraviolet logarithmic behavior. Also note that earlier truncations 
of the ghost-gluon system, with kernels violating multiplicative renormalizabil- 
ity, yield a very different infrared analysis where only the ghost loop contributes 
to the leading infrared behavior of the gluon vacuum polarization. 

The infrared analysis of the one-loop diagrams will be done in the same 
way as before (see ref. [25] and references therein), and the method will be 
extended to compute the infrared contributions of the two-loop diagrams. As 
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will be detailed further on, both gluon and ghost equations are individually 
satisfied by propagators obeying infrared power laws 

i^K(x) = Ax«, Gr{x) = Bx^, (50) 

where x = p^. When substituting these power laws in Eqs. (19) and (20), the left 
and right hand sides of the integral equations yield power behaviors. Requiring 
the leading infrared exponents on both sides to match, produces the constraint 
a = —2/3, such that 

Fr{x) = Ax''^, Gr{x) = Bx-''. (51) 

Because of Eq. (35), these power laws generate an infrared fixed point for the 
running coupling 

ao = lim aij,FR{x)GR{x) = a^AB^ constant . (52) 

At this point it is useful to introduce the related function 

Hx) = = -^Fr{x)Gr{x) . (53) 

Prom Eqs. (19), (20), (24), (42)-(46) it is clear that the function \{x) only de- 
pends on its value at the renormalization point, and is independent of the 
number of colors N^. for the MR truncation presented herein^. As will be shown in 
this section, the infrared value Ao is fixed by the coupled Dyson-Schwinger equa- 
tions. Moreover, as Aq is independent of iVc, Eq. (53) shows that the infrared 
fixed point ckq of the running coupling is proportional to 1/Nc- 

Equating the coefficients of the leading power of x, for x ^ 0, on both sides 
of each equation yields two relations between the infrared fixed point Aq and the 
leading infrared exponent k: 

/g(Ao,k) = , /ir(Ao,K)=0. (54) 

The existence of a consistent infrared power solution requires the gluon and ghost 
equation conditions (54) to be satisfied simultaneously. The solutions of Eq. (54) 
will determine the values of the infrared fixed point Aq and exponent k. The 
various loop integrals will now be computed and discussed. 

4.1 Ghost Equation 

The leading power contributions to the ghost equation are calculated by substi- 
tuting^ the power laws (51) in Eq. (19) with self-energy (42): 

^ = Zs- N^glAB ^ To(x, y, z) y^^ z^^ , (55) 

*Note that this is also true up to third order in perturbation theory [32]. At four- loop order 1/A^c 
corrections are generated [33] due to non-vanishing contributions from non-planar diagrams. 

®Even though the general solutions behave logarithmically at high momentum {x ^ Aqcd'^), 
power laws are used for all momenta in the infrared analysis. This procedure gives the correct 
leading power contribution because pure power laws are themselves solutions of the coupled 
equations in the MR truncation, corresponding to Aqcd — > oo. Ref. [13] showed that the 
leading infrared power is independent of Aqcd, which only enters the solution through the 
next term of its infrared asymptotic expansion. 
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where x = p^, y = Q^, z = {p — q)"^, k = {ki, ^2) is defined by k = K(2b£ — a^), 
and ki + k2 = K because as,i + 0^,2 = &s,i + ^e,2 = 1 (from Eqs. (41) and (42)). 
Substituting the kernel Tq given in (47) yields 



= BZ3 + NcglAB^ djm J x' y 

i.i.m. ^ 



i,i,m 



where i.j, m, arc integers satisfying i + j + m = —2, and the coefficients 

be read off from Eq. (47). The right hand side is a sum of integrals of the form 



x'^yl'z'^, (57) 

with a + /3 + 7 = k — 2. In ref. [20] these integrals were solved using spherical co- 
ordinates, and the results were expressed in terms of generalized hypergeometric 
functions. A more concise, equivalent expression can be derived [22] by noting 
that integrals of type (57) are typical scalar integrals encountered when apply- 
ing dimensional rcgularization in perturbative calculations and they are readily 
computed using Eq. (D.5). 

Equating the coefficients of the leading infrared power x'^ on both sides of 
Eq. (56) gives 

XoXG{ki,k2) = l, (58) 

with Ao = NcaQ/4:TT and cto given by Eq. (52). The function XG is given by the 
self-energy integral in Eq. (56) which can be calculated analytically using Eq. 
(D.5). Its general expression is 

, . 3 r(2 + fci)r(l + fc2)r(-A:i-A:2) 

= -2r(i-W2-/c2)r(3 + fci + M ' 

where the Gamma function recurrence relation was used repeatedly to bring the 
expression into its simplest form. The values of ki and /c2 depend on the details 
of the truncation through the sets as and bs of (42) which have not yet been 
fixed. 

For the ghost self-energy no symmetry between internal legs can be used to 

fix as and bs, and three different truncations are considered. A natural choice 

corresponding to the bare ghost-gluon vertex approximation is obtained with 

IR 

as = (1,0) and bs = (0,1), yielding a kernel proportional to GR{y)Fji{z) ^ 
y-K^2K_ rpjjjg (.];ioice glvcs {ki,k2) = {—k,2k), and from Eq. (59) 

n / X / X 6r2(2-K)r(2K) 

An expression equivalent to Eq. (60) was first given in ref. [20], where it was 
written in terms of generalized hypergeometric functions. 

The G/G nature of the ghost-gluon vertex corrections makes it possible for 
the bare vertex to satisfy MR and reproduce PT. Nevertheless, it is clear that 
the exact dressed ghost-gluon vertex is not identical to the bare one, as can for 
example be seen from its Slavnov-Taylor identity [12]. That identity suggests 
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vertex corrections like G{r'^)/G{q^) where q is the ghost momentum in the loop, 
and r the gluon momentum. Approximating the full ghost-gluon vertex by a bare 
vertex multiplied by such a correction leads to the truncation as = (0, 1) and 
bs = (0, 1), and the kernel is proportional to GR{z)Fji{z) with infrared limit z'^, 
i.e. ki = 0, k2 = K- In this case 

XU-) - Xg{0, .) = 2,(1 + ,)(2 + • (61) 

However, as explained in the previous section, in general the MR truncation 
requires an additional vertex type correction to satisfy MR and reproduce the 
correct perturbative behavior of the integrals. This is not obvious for the ghost 
self-energy because of the G/G nature of the vertex corrections, but is clear for 
diagrams containing three-gluon, four-gluon and quark-gluon vertices. Introduc- 
ing such an additional correction leads to another variant, which is very similar 
to the multiplicatively renormalizablc truncation of the quark equation [34] . This 
corresponds to as = (—1,2) and bs = (0, 1), yielding a kernel proportional to 
a{z)/Gji{y) (using Eq. (35)) with infrared limit y'^, i.e. ki = k, k2 = 0, and 

xSW.X«(-<,0) = ^;^. (62) 

The three truncations (60), (61) and (62) will respectively be referred to as 
bare truncation, V-truncation and a-truncation. The variation of xg function 
of K for these three cases is shown in Fig. 1. The plotted quantity is inversely 
proportional to the infrared fixed point of the running coupling (see Eq. (58)), 
and from the figure one can get information about the dependence of the fixed 
point on k using different kernel corrections. For instance, a value of k = 0.5 
would correspond to a fixed point value Aq = 5/8 for the bare truncation and 
V-truncation, and Aq = 5/6 for the a-truncation. Global infrared and ultraviolet 
convergence criteria of the integrals require < k < 2. 



4.2 Gluon Equation 

Next, the various contributions to the gluon vacuum polarization are investi- 
gated. It is important to note that, in contrast to the studies [12, 13, 20], all the 
vacuum polarization diagrams (43) to (46) of the MR truncation will be shown 
to contribute to the leading order infrared power of the gluon equation (20). The 
infrared analysis yields the following equation after substituting the power laws 
(51) in Eq. (20) 

^ = Z3- lim[n2'=(x) + n2a(a;) + U^a^x) + n^«(x)] , (63) 

or for the coefficients of the leading power in x: 

1 = - lim Ax'^'' \U'^%x) + U^^ix) + n3»(x) + n4»(x)l . (64) 
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0.2 0.3 0.4 0.5 0.6 0.7 0.8 

K 

Figure 1. Variation of the infrared ghost self-energy integral xc as function of the infrared 
exponent k for three specific choices of the kernel in Eq. (42): Xg (60) using a bare ghost-gluon 
vertex, Xa (61) with a factorizable vertex Ansatz, and Xg (62) for a a{z)/GR{y) kernel. 



Each contribution is now investigated in turn. The infrared ghost loop con- 
tribution (43) to Eq. (64) is given by 



t,j,m - 1 



(65) 

with X = p'^, y = q'^, z = {p — q)^. The integers i,j,m satisfy i + j + m = —2, 
and the coefficients Cijm can be read off from Eq. (48). Furthermore, k = (ki, ^2) 
is defined by k = K(2b2c — ^2c), and ki + k2 = —2k because a2c,i + a2c,2 = 2, 
&2c,i+?'2c,2 = (fromEq. (43)). 

Although one can easily calcTilate this contribution for arbitrary a2c and b2c, 
their values are naturally fixed by symmetry considerations between the two 
internal ghost legs: a2c = (1, 1), b2c = (0,0). The kernel is then proportional to 
GR{y)Gji{z), and its infrared contribution, for which k = (— k, —k), is 

, _ r^(2-K)r(2K) 

This choice corresponds to an effective bare ghost-gluon vertex approximation. 
Although one could argue that different ghost-gluon vertex approximations are 
used in the ghost (Sect. 4.1) and gluon equations, this is not an inconsistency 
as the truncation considers effective kernel corrections rather than Ansatze for 
the full dressed vertices. Because of the different tensor structures of the ghost 
self-energy and gluon vacuum polarization diagrams, kernel corrections can be of 
different natures for both equations. Furthermore, for the ghost loop contribution 
to the gluon vacuum polarization, corrections analogous to those of Sect. 4.1 
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would depend on the external gluon momentum, and are not allowed in the 
effective kernel truncations developed in this paper as they would destroy the 
perturbative properties of the truncation. 

The evolution of is shown in Fig. 2 together with the ghost curves for 
Xq, Xg Xg- the leading ghost-loop truncation of ref. [20] the consistency 
condition (54) between ghost and gluon equations is I/Aq = Xg{k) = Xf ('^)) a-^d 
the crossing of x|f and Xg at k = 1, with Aq = 1, was suggested as a possible 
solution of the coupled equations. Ref. [20] only described the infrared analysis, 
and the numerical study of the equations did not reveal solutions satisfying 
the equations over the whole momentum range. Using the a-truncation (62) 
of the ghost equation would yield a slightly different infrared solution in the 
leading ghost-loop approximation, with k = 0.79 and a relatively large coupling 
Ao = 1.46. Note that the results can be sensitive to the truncation as can be 
seen from the V-truncation for which no solution for k exists, as no intersection 
between this ghost curve and the gluon curve occurs in Fig. 2. 



Figure 2. Variation of tlie infrared 
ghost loop integral Xf (66) in the gluon 
vacuum polarization, together with the 
ghost solf-oncrgy integrals Xg (60), Xg 
(61) and Xg (62) from Fig. 1, as func- 
tion of the infrared exponent k. The 
crossing of Xf^ and Xg at k = 1 was 
suggested as a possible solution for 
the equations in the leading ghost-loop 
truncation of ref. [20]. 

The full numerical investigation of the leading-ghost loop truncation was 
repeated in the current study and it confirms that no solution can be found 
connecting the infrared to the ultraviolet asymptotic behavior. The failure to 
find numerical solutions does not prove that such do not exist, as the numerical 
solution of the equations could be sensitive to instabilities. However, as will be 
shown later on, the numerical method applied to solve the equations is quite 
powerful, and it is therefore very likely that this truncation has no solution at 
all, and that the infrared solutions are spurious solutions of the infrared analysis. 
This problem will be resolved in an elegant manner by the MR truncation. 

In a further attempt to find numerical solutions of the coupled equations in 
the leading ghost-loop approximation with exact angular integrals, refs. [22, 23, 
24] extracted the gluon dressing function from the tensor equation (4) using V^v 
of Eq. (22) with arbitrary n. They reported new infrared solutions with k ^ 0.5 
for n 7^ 4. These solutions are spurious as is explained below. 

Although the solution can depend on n if the DSE truncation violates the 
transversallity of the gluon vacuum polarization, one hopes that its sensitivity 
is minimal. However, an additional problem occurs when using a non-BP tensor 
(n 7^ 4) as spurious quadratic divergences are generated in the gluon vacuum 
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polarization [10]. These divergences have to be removed properly in order to get 
renormalizable solutions. Performing the infrared analysis of the ghost loop in 
the gluon vacuum polarization for arbitrary n leads to an extension of Eq. (66) 

2e,„ ..^_r„ , o , , . , r2(2 - k) r(-l + 2 «:) 



XF{n,K) = [n (-3 + 4k) + 6-4k] 



6r(4-2K) r2(l + K) 



= ,|»(4,„)+^(-3 + 4K)g|^^fIti+M, (67) 

where x|^^(4, k) is given by Eq. (66). Equation (67) shows how ^^^{n^n) can be 
written as a sum of two terms: one regular and the other singular at k = 0.5. 
It is easy to check that the singular factor r(— 1 + 2k) is exactly generated by 
the terms of Mq (extension of Eq. (48) for arbitrary n) which are quadratically 
divergent in the ultraviolet limit. The singularity only occurs for n 7^ 4 and is 
spurious. Figure 3 shows x|f for n = 4 and n = 3, together with the ghost curve 
Xq. The singularity at k = 0.5 for n 7^ 4 inevitably introduces an additional, 
spurious crossing of the ghost and gluon curves just above k = 0.5. 



4 




I ' ^ ' ' ' 1 

0.2 0.4 0.6 0.8 1 1.2 

K 

Figure 3. Variation of the infrared ghost loop integral xl? (67) in the gluon vacuum polarization 
for n — A and n = 3, together with the ghost self-energy integral Xc; (60), as function of the 
infrared exponent k. The additional crossing for the n = 3 curve is spurious as it is due to the 
divergence of Xf^ at k = 0.5, which is related to the presence of spurious ultraviolet quadratic 
divergences when n ^ 4 in Eq. (22). 

When removing the kernel terms responsible for the quadratic divergence 
in a natural way, the singular Gamma function in Eq. (67) disappears and no 
spurious solution for k is found. The only solution left for n 7^ 4 is the one that 
continuously connects to the k = 1 solution for n = 4 [20]. Unfortunately, after 
removal of the quadratically divergent terms the integral equations do not have 
overall solutions, as was the case for n = 4 in ref. [20], and as could be expected 
the absence of solutions is independent of n. 
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The absence of numerical solutions is remedied in ref. [24] by introducing 
a specific procedure to remove the quadratic divergences for n 7^ 4. The terms 
responsible for the quadratic divergences are not simply removed from the respec- 
tive ghost and gluon loop kernels where they originally occur. Rather, the sub- 
traction in the gluon loop is adapted so that it cancels the ultraviolet quadratic 
divergences coming from both gluon and ghost loops, while the ghost loop is left 
unsubtracted. As the gluon loop is infrared subleading in that truncation, the 
subtraction scheme leaves the leading infrared analysis of Eq. (67) unchanged. 
The truncation yields infrared solutions with k ^ 0.5 for n 7^ 4, as indicated in 
Fig. 3, and is free of ultraviolet quadratic divergences at the same time. Further- 
more, the numerical treatment found overall solutions for the integral equations, 
connecting the asymptotic behaviors at low and high momenta. 

This truncation scheme was the first one to yield full numerical solutions in 
a treatment with exact angular integrals (for any n, except n = 4). Nevertheless, 
it is difficult to evaluate the usefulness of the truncation because of the pecu- 
liar scheme used to subtract the quadratic divergences. Clearly, this subtraction 
scheme can only be used as long as the gluon loop does not contribute to leading 
infrared order, and this is in contradiction with the MR truncation of Sect. 3. 
The goal herein is to go beyond that truncation, avoiding quadratic divergences 
from the very beginning by using n = 4, and to develop a truncation respecting 
multiplicatively renormalizability, for which all diagrams contribute to leading 
infrared order. 

As was shown in ref. [25] the gluon loop also contributes to the leading 

infrared power behavior of the gluon vacuum polarization in the MR truncation. 
The treatment is analogous to that of the ghost loop, and its contribution is 
given by 

Ao xl'{ki,k2) = - lim ^x^^tf ^(x) , (68) 

with k = {ki,k2) defined by k = K(2b2p — &2g), and fci + /c2 = —2k because 
02g,i + ^25,2 = 2, 62^,1 + b2g,2 = 0. This Icads to a sum of integrals similar to that 
of Eq. (65), where the coefficients Cijm are now determined by the kernel (49). 

The values a2g and h2g are again naturally fixed based on symmetry con- 
siderations between the two internal gluon legs: a.2g = (1,1), h2g = (0,0). The 
kernel is proportional to GR{y)Gii{z), and its infrared contribution is^ 

Note that in this case the effective MR kernel automatically corresponds to some 
implicitly defined dressed vertex. Using a bare triple gluon vertex in Q of Eq. (38) 

would yield equations that do not satisfy Eq. (44), would violate multiplicative 
renormalizability and would furthermore not contribute to the leading infrared 
power. 

When neglecting the two-loop diagrams and only taking into account the 
one-loop contributions in the gluon vacuum polarization, the infrared constraint 

®The expressions (66) and (69) given here in terms of Gamma functions are identical to 
those expressed in terms of generalized hypergeometric functions in ref. [25], as was checked 
numerically. 



J.C.R. Bloch 



19 



(64) from the gluon equation is 

Ao[xlf(«)+XF'W] = l- (70) 

The existence of consistent infrared power solutions requires the ghost and gluon 
relations (58) and (70) to be satisfied simultaneously: 

XG{^) = xfi^)+xl'{'^)- (71) 



Solutions of this equation would yield possible values for the leading infrared 
exponent k. The corresponding ag would then be computed from Eq. (58) or 
(70). In Fig. 4 the gluon loop contribution Xp is drawn together with Xf 
Xg, and it is clear that the condition (71) cannot be satisfied for any k between 
and 1 because of the large negative contribution of the gluon loop. This lead 
to the conclusion in rcf. [25] that the MR one-loop-only truncation, where only 
the ghost and gluon loops are retained in the gluon vacuum polarization, does 
not support infrared power laws for the propagators. 




Figure 4. Contributions Xg: xP S'Hd 
X^/ as function of the infrared expo- 
nent K. The ghost curve Xg and the 
gluon curve Xf+Xf '^'^ ^'^^ intersect, 
and power laws are not solutions in 
the one-loop-only MR truncation. 



The absence of power law solutions in that truncation could be interpreted as 
a setback, but it was suggested in ref. [25] that taking into account the two-loop 

diagrams in the DSEs could substantially change the infrared behavior of the 
equations. Herein it will be shown that the absence of consistent infrared power 
law solutions is an artifact of an incomplete truncation due to the omission of 
the two-loop diagrams, and that taking these into account, albeit in a simple 
realization, will reinstate power laws in a very natural way. 

The following section will illustrate how the infrared analysis of the two-loop 
gluon vacuum polarization diagrams can be performed quite elegantly. 



4.3 Two-loop Diagrams 

The sunset contribution (45) to the gluon vacuum polarization in the MR trun- 
cation is given by 

Ii^^{x) = Nlglj [ Q0(p,gi,g2,g3)G|^[a3,]FA[b3,]. (72) 
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Taking the infrared limit yields the contribution of the sunset diagram to the 
coefficient of the leading infrared power x~^'^ in the gluon equation (64): 



xl'ih, k2M) = - lim ^x2-n33(a 



(73) 



Substituting the power laws (51) yields^ 

Xf'(^i> k2M) = -(167r2)2 / / Ql(j>, q^,q2, q^) x^'^ z'^^ z^^ , (74) 

./ 01 Jar> 



with X 



qf. Furthermore, k = (fci, k2, k^) is defined by k = K{2h 
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asg), and ki+k2+k3 = -2k because 033,1+033,2+033,3 = 4, bsg^i+b3g^2+hg,3 = 1- 
To compute the sunset diagram it is crucial to note that it can be considered 
as a sequence of two single loop integrations. First, the (g2, '73)-loop, enclosed in 
the box in the diagram shown below, is integrated out, and then the remaining 
{qi,qs)-loop, with qs = q2 + = P - Qi, computed. 







92 























(75) 



93 



The integration rule (D.5) cannot be applied as such to compute the inner loop, 
because the kernel Q3 cannot be written as a function of squared momenta only, 
as could be done for the one-loop diagrams. There is always one angle that cannot 
be eliminated. Fortunately, the integral can still be solved in a fairly easy way 
by treating it as a sum of tensor integrals of rank 4 or lower. 

The various tensor terms in the kernel of Eq. (74) have denominators pro- 
portional to (h^^q^^Q^^i with non-integer powers 01,02,03. To transform this 
into known integrals, first the Feynman trick (C.l) is applied to rewrite the 
denominators of the inner loop integral as 



1 



r(02 + 03) t«2-l(l_i)a3-l 



ql'^'qf' r(02)r(03) Jo [tql + (1 - t)qf[ 



(76) 



The momentum q^ is eliminated by introducing the outer loop momentum qg 

92 + 93 



1 



r(a2 + 03) 



dt 



ql'^^qf- r(02)r(03) 70 [ql - 2(1 - t)q2 ■ + (1 - t)?|] ' 



(77) 



To eliminate the dot product in the denominator, so that it is written in terms of 
squared momenta only, the variable transformation £ = q2—{l—t)qs is introduced: 



r(o2 + 03) 



ql^'ql"' r(02)r (03) Jo [^2 + _ t)g2]«2+a3 " 

'^Power laws are assumed for all momenta as explained in footnote 5 on page 12. 



(78) 
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The inverse variable transformations 

q2=e+il-t)qs , q3=tqs-e, (79) 

are then used to replace the four-momenta components of q2 and qs by i and 
qs in the numerators of the tensor integrals of Eq. (74). The ensuing inner loop 
integrals can be calculated using the tensor integral formulae Eqs. (D.l) to (D.4). 
Consequently, the remaining integration over the Feynman parameter t is easily 
performed using Eq. (C.3). After integration of the inner loop the kernel of the 
external {qi, qs) loop integral is a function of p, qi and qg which can be written 
as a function of squared momenta only, and this integral can be computed using 
Eq. (D.5). This complete procedure was programmed using the symbolic algebra 
program Form [35] and yields the following result for the integral (74) with Qg 
defined as the tree-level version of kernel (33): 

Xl^iki,k2,k3) = l{l + ki+k2 + k3-2kik2-2kik3-2k2k3-8kik2 fes) 



r(i + k^) r(i + k2) r(i + k^) rj-k^ -k2- fcg) 
r(2 - ki) r(2 - k2) r(2 - ^3) r(4 + ki + k2 + k3 



(80) 



As expected from diagram (75), this result is invariant under ki ^ k2 ^ /C3. The 
parameters in the MR truncation of the sunset diagram (72) are chosen to respect 
the symmetry between the three internal gluon legs, i.e. asg = (4/3, 4/3, 4/3) and 
hsg = (1/3, 1/3, 1/3) corresponding to a kernel 

Qlip, 91, 92, 93) Fl,/\z,)Gf{z,)F'J\z2)Gfiz2)Fi/\z3)Gf{z3) . (81) 

Hence, the set k in the infrared equation (74) respects this symmetry too, and 
ki = k2 = k3 = —2k/3. For such a choice, Eq. (80) becomes 

X%'{k) = X^^(-2«:/3, -2k/3, -2k/3) 

80 64 o\ r3(i-^)r(2K) , , 

1-2k--k^ + —k^ , ^ o ^ . \ ■ (82) 
3 27 y r3(2 + f )r(4 -2k) ^ ^ 

The treatment of the squint diagram can be performed in a way analogous 
to that of the sunset diagram. The squint contribution (46) to the gluon vacuum 
polarization in the MR truncation is 

n^^(x) = N^gf, [ [ QHp, qi,q2, 93, 94) F^K] . (83) 

Jqi Jq2 

Again, the contribution to the coefficient of the infrared leading power is 
given by taking the infrared limit, 

X%'{ki, k2, k3, ki) = - hm Ax2«n^s(x) , (84) 
and substituting the infrared power laws: 
XF'{ki,k2, k3, ki) = -(167r2)2 [ j Ql{p, q,,q^, q,, q^) x^'' z\^ z''^ zf , (85) 
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with X = p^, Zi = qf. Furthermore, k = {k\,k2, k^, k^) is defined by k = K{2h,ig — 
a4g), and ki + k2 + k^ + k4 = -2k because 043,1 + 043,2 + aAg,3 + 043,4 = 4, 

As before, the two-loop integral can be computed by evaluating two consec- 
utive one-loop integrals: first the inner-loop integral over {q2,q3) and then the 
remaining outer- loop integral over (gi, 54) as shown in the diagram below. 




The steps to be followed are exactly the same as those described above for the 
sunset diagram. The inner loop integral over (^2,^3) can be transformed into 
known tensor integrals of rank 5 or lower. Once these have been integrated out, 
the remaining one- loop integral has a kernel which is a function of p, qi , q^ which 
can easily be solved using Eq. (D.5). The final result is 

Xf (^1)^2,^3,^4) = -- (3/c2^ + k2 k-i + Sks'^ + k2 + k3 - 5) 

o 

X (13 + 3 (fei + k2 + ks) + 13 ^4 + 6 (fci + fcs + ^3) ^4 - 6 ^4^) 

^ r(i + fc2)r(i + fc3)r(-fc2-A:3) 
r(2 - k2) r(2 - ks) r(3 + k2 + k^) (87) 
r(i + kj) r(2 + fci + fca + fcg) r(-fci - fc2 - fcg - fc4) 

^ r(2 - k^) r(l -ki-k2- ks) r(4 + ki + k2 + k3 + ki) ' 

Note that is symmetric under interchange of q2 and gs as is expected from 
diagram (86) . The value of Xj? depends on the choice of a.4g and in Eq. (46) . 

In the next section, the existence of infrared power solutions will be investi- 
gated taking into account all the contributions calculated above. 



4.4 Consistency Condition 

The infrared consistency condition (71) has to be modified when the two-loop 
diagrams are taken into account in the gluon vacuum polarization. For the gluon 
equation the relation between the infrared fixed point of the coupling and the 
infrared exponent of the propagator power laws can be derived from Eq. (64) 
using Eqs. (65), (68), (73) and (84): 

l = AoXFW + AgxF(«), (88) 
where the one-loop and two-loop contributions are defined by^ 

x¥{^)=x¥{^)+xl'{^) , Xfi'^) = xl'i'^) + XF'{^)- (89) 
*The x's depend on k = K{2b — a), and can be written as functions of k for fixed a's and b's. 
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Equation (88) is solved by 

^ = \ {x'A^) ± \l[x'Mf + ^x'^{^)) . (90) 

As discussed in Sect. 4.2, the one-loop contribution is negative, and therefore a 
positive solution for Aq can only occur for the "+" -solution in Eq. (90). 

The infrared consistency condition between ghost and gluon equations simply 
follows by requiring Eqs. (58) and (90) to be satisfied simultaneously, and the 
infrared exponent k, is determined from 

Xg(«) = \ (x'A^) + V[x^'('^)]' + 4xf('^)) ■ (91) 

For the simple MR kernel truncations developed in this paper the 
parametrization of the ghost loop (43) , gluon loop (44) and sunset diagram (45) 
were fixed by symmetry considerations between the internal momenta. However, 
the parametrization of the squint diagram (46) has not yet been fixed as a4g, 
still have to be assigned specific values. A major concern of this study is to in- 
vestigate if sets can be found for which the infrared consistency condition (91) 
is satisfied, and if so, if numerical solutions for the propagators exist having this 
asymptotic infrared behavior. 

First, the situation is investigated where all diagrams are taken into account 
except the squint diagram. The evolution of X/? (82) as function of k is shown 
in the left pane of Fig. 5, and for example x^(0.5) ~ —0.1. Table 1 summarizes 
the various infrared contributions in Eq. (91), with exception of x^^, for k = 
0.25,0.5,0.75, and in the right pane of Fig. 5 the evolution of 1/Aq from the 
gluon equation condition (90) is shown. The *-superscript is used to emphasize 
that the squint diagram was neglected. Clearly, no intersection will exist between 
the gluon curve of Fig. 5 and any of the ghost curves of Fig. 1, and no power 
solution exists which simultaneously satisfies both ghost and gluon equations. 

Within the scope of the MR truncation, the only possible way left to find a 
consistent power solution to both equations is for the squint diagram to give the 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 



K K 

Figure 5. Left pane: variation of the sunset contribution Xf* (§2) as function of the infrared 

exponent k. Right pane: Inverse of the infrared fixed point \q (90) as function of the infrared 
exponent n without squint contribution. 



24 



Two-loop Improved Truncation of the Ghost-Gluon DS Equations 



Table 1. Infrared contributions from the gluon vac- 
uum polarization diagrams and ghost self-energy. An 
infrared solution requires 1/AS from the gluon equa- 
tion to equate xg from the ghost equation. 



K = 0.25 K = 0.5 K = 0.75 



gluon Xf(k) 


0.55 


0.5 


0.65 




-6.6 


-5.4 


-7.3 




0.08 


-0.1 


-0.60 


1/AS(k) 


0.01 


-0.02 


-0.09 


ghost xg(«;) 


2.92 


1.6 


1.21 


Xg(«) 


2.67 


1.2 


0.73 



contribution necessary to balance both sides of Eq. (91). A consistent infrared 
solution requires the squint diagram contribution to satisfy 



?9l 



Jig I 



(92) 



and its variation with k is shown in the left pane of Fig. 6. The figure shows that 
the squint contribTition is quite constrained for a solution to exist: it has to be 
positive and larger than 5. The figure also illustrates that the sunset diagram is 
negligible in the analysis of the infrared asymptotic behavior. 

As discussed in Sect. 3, the MR truncation assumes the cancellation of a cer- 
tain number of dressing function corrections in the loop kernels. The parameters 
a4g, hiig of Eq. (46) determine how these corrections are distributed over the 
internal legs of the squint diagram. Five simple choices are shown in Table 2 as 
illustrations. Set 1 is a simple equipartition of the dressing corrections over the 
various legs. Sets 2-5 correspond to attaching a G/i^-correction to each three- 
gluon vertex and G"^ /F to the four-gluon vertex. The sets differ in the way the 
dressing corrections are distributed over the vertex legs. 




50 





/ ' /set1 L 




/ / set 2 i- 




/ sets ; 


\ \ 


/ set 4 \ 




I sets \ 



0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
K 



Figure 6. Left pane; Required value Xf* (92) as function of the infrared exponent k to ensure 
consistent infrared power solutions: (A) using the bare ghost equation truncation Xg (60) , (B) 
for the a-truncation Xg (62). The dashed curves, calculated without the sunset contribution, 
show that its contribution is negligible. Right pane: Squint contribution Xj? (87) as function of 
the infrared exponent k for the parametrizations of Table 2. 
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Table 2. Various truncations of the squint diagram. 



st^t 



a.i,, 



b4„ 



k = K(2b_i,, - a.i.,,) 



1 (1,1,1,1) 

2 (5/6,1,1,7/6) 

3 (5/6,4/3,4/3,1/2) 

4 (1/2, 7/6, 7/6, 7/6) 

5 (1/2,3/2,3/2,1/2) 



(1/4, 1/4, 1/4, 1/4) 
(1/6,1/3,1/3,1/6) 
(1/6,1/6,1/6,1/2) 
(1/2,1/6,1/6,1/6) 
(1/2,0,0,1/2) 



(-k/2, -k/2, -k/2, -k/2) 
(— k/2, — k/3, —k/3, —5k/6) 
(— k/2, — k, —k, k/2) 
(k/2, -5K/6, -5«:/6, -5k/6) 
(k/2, -3k/2, -3k/2, k/2) 



As can be seen from Eq. (87) the infrared contribution of the squint diagram 
depends on k = K(2b4g — a4g), rather than on a4g,b4g themselves, and the 
fc-partitionings corresponding to the various parametrizations are also given in 
Table 2. The infrared contributions for these squint parametrizations are 
shown in the right pane of Fig. 6. Comparing the required and actual squint 
contributions in the left and right panes of Fig. 6 demonstrates that the squint 
diagram potentially has the right sign and magnitude to satisfy the infrared 
consistency conditions. Therefore, the value of k and Aq can be determined by 
looking for the intersection between the ghost and gluon curves. This is illustrated 
in Fig. 7 for one specific squint truncation, i.e. set 4 of Table 2. Mathematica 
was used to determine the solutions of various truncations of the ghost-gluon 
system, constructed by combining all the squint truncations of Table 2 with 
two different truncations of the ghost self-energy. Each truncation has up to 
two infrared solutions which are tabulated in Table 3. The numerical analysis of 



^ 3 



Gluon eq. - set 4 
Ghost eq. a-trunc. 




0.1 
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0.8 
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Figure 7. Inverse fixed point 1/Ao from the gluon equation, with squint parametrization set 4 

of Tabic 2, and from the ghost equation in a-truncation, as function of the infrared exponent k. 
The intersection between gluon and ghost curves gives the solutions of the infrared consistency 
conditions. Two solutions are found: k = 0.36 and 0.96. Only the first one corresponds to a 
true solution of the integral equations, as was established from the numerical study. 
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Table 3. Solutions of the infrared consistency condition (91) and exis- 
tence of numerical solutions for various truncations of the ghost-gluon 
system. Part A for the bare ghost self-energy truncation and part B for 
the a- truncation. The squint truncations 1-5 are defined in Table 2. 

truncation solution 1 solution 2 



ghost 


gluon 


K 


Ao 


numerical 




Ao 


numerical 


A 


1 


0.24 


0.32 


yes 


0.95 


0.96 


no 




2 








0.78 


0.84 


no 




3 


0.38 


0.50 


yes 










4 


0.53 


0.65 


yes 


0.90 


0.92 


no 




5 


0.35 


0.46 


yes 


0.90 


0.93 


no 


B 


1 


0.17 


0.25 


yes 


0.98 


1.95 


no 




2 








0.91 


1.77 


no 




3 


0.31 


0.48 


yes 










4 


0.36 


0.57 


yes 


0.96 


1.89 


no 




5 


0.31 


0.48 


yes 


0.93 


1.81 


no 



the coupled integral equations, which will be discussed in more detail later on, 
showed that not all infrared solutions correspond to consistent overall solutions 
connecting the infrared and the perturbative ultraviolet behaviors. The solutions 
with K ^ 0.75 are spurious solutions of the infrared analysis, and only those with 
K G [0.17,0.53] correspond to consistent solutions of the integral equations. 

Clearly, the simple MR truncation is not able to pin down the exact QCD 
result: the value of n depends on the truncation and its range is quite broad. 
Nevertheless, the study elucidates the conspiracy between the various diagrams 
in the propagator equations, necessary to allow for consistent solutions of the 
coupled ghost-gluon system. 

Now that power law solutions have been found satisfying the infrared consis- 
tency conditions, the next task is to solve the integral equations numerically and 
look for overall propagator solutions having such an infrared behavior. However, 
first a closer look is taken at the numerical treatment of the two-loop diagrams. 

5 Numerical Study 

5.1 Approximated Two-loop Diagrams 

Although the infrared analysis of the two-loop diagrams could be performed ana- 
lytically, their numerical treatment as part of an integral equation is a numerical 
problem of too high complexity to be accessed yet. Therefore an alternative so- 
lution method is herein proposed where the two-loop integrals arc replaced by 
effective one-loop integrals which have the same leading infrared behavior as the 
original integrals, are subleading in the perturbative region, and still satisfy the 
properties of multiplicative renormalizability. 

The MR truncated squint diagram was given in Eq. (46) by 

U^^ix) = Nlgl [ [ Qlip, qi,q2, qs, q,) G|j[a4,] F^[b4,] , (93) 
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where a4g = (ai, 02, 03, 04) and h4g = (61, 62, ^3,^4) parametrize the squint trun- 
cation. 

As discussed in the infrared analysis (see Sect. 4.3), the two-loop integration 
can be performed in two successive steps: 

U'^{x) = N^gf, [ Utljx,z,,z,) G''j^{z,)f'^\z^)G%\z,)f'^\z,), (94) 

Jqi 

with inner loop integral 

Ii^Ux,zi,z,)= [ Ql{p,qi,q2,q3,q4)G'g{z2)F'^\z2)G';i{z3)F'^%Z3). (95) 

This two-step procedure provides an ideal opportunity to introduce effective 
one-loop integrals as approximation to the two-loop integrals, based on their 
analytic infrared evaluation. To evaluate the infrared contribution of the squint 
diagram, the dressing functions are replaced by their asymptotic power behav- 
iors, and the inner loop (95) can be written as 

Jq2 

with k2 = (262 — 0,2)^ and ks = (263 — 03)^. The integrals are evaluated as 
detailed in Sect. 4.3, and the solutions can be written as a function of squared 
external and outer loop momenta: 

lim nf„^^^^(x, Zi,Z^) = ^f2+635a2+a3^fe+fc3 ^ ^^^.^(^^^ ^ (97) 

i,j,m 

where the coefficients Cijm are ratios of F- functions depending on A;2, ^3. 

The powers in front of the sum can be identified with the infrared behaviors 
of Gii{zi) and Fii{zi) and, in the infrared, the last expression is identical to 

1™ nfn'ner(^> ^1, ^4) = (^1)^^^+^^ (^1) C,jm{k2,ks) 4zT • (98) 

i,j,m 

Although Eq. (98) is only exact in the infrared limit, it will be used to ap- 
proximate the inner-loop integral (95) for arbitrary external momentum x, in 
order to construct an effective one-loop diagram replacing the original two-loop 
diagram (94) in the numerical calculations. Substituting Eq. (98) in Eq. (94) 
yields the following effective squint contribution: 

n^^(x) ^ iV^^ / Qf{x, y, z) G"^ (y) f'^ {y) G^^ (z) F'^ {z) , (99) 

Jq 

with 

= ai -f 02 + as , = 04 , ^^^^^ 

hy = hi + h2 + h, bz = b4, 

and where the effective kernel 

Qf{x, y,z)=J2 Cijm{k2M) x'y'z"" , (101) 
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for arbitrary external momentum x, is defined by the infrared integral of the 
inner loop, given in Eq. (97). Its explicit form is given in Appendix E. 

Because of the properties of a.4g and h4g, the exponents in the effective loop 
obey 

ay + a^ = A by + b^ = l, (102) 

so that the approximation conserves the number of dressing functions of the 
original integral. This is sufficient to ensure that the approximation remains 
multiplicatively renormalizable and ultraviolet subleading. Furthermore, the ex- 
ponents also satisfy 

(2by - ay)n = ki + k2 + kz , {2bz — az)K = , (103) 

and the effective loop reproduces the same infrared result (87) as the original two- 
loop diagram. The preservation of these properties is obvious by construction. 

Away from the infrared limit, the propagators no longer behave as power laws, 
and therefore the infrared powers in Eq. (97) were replaced by the corresponding 
dressing functions in Eq. (98) when devising the approximation. This guarantees 
that the ultraviolet behavior is logarithmically subleading and avoids contamina- 
tion with unjustifiable power behaviors. Indeed, substituting the leading-order 
perturbative results F{x) ~ a{x)'''^ and G{x) ^ a{x)^^ , with 7 = —13/22, 
5 = —9/44, in the approximated two- loop diagram (99) yields an ultraviolet con- 
tribution of order a(.x)^^'^, which has the correct log-exponent of the resummed 
next-to-lcading order correction to the gluon anomalous dimension [36]. Hence, 
the effective loop correctly mimics the ultraviolet behavior of the two-loop dia- 
gram, albeit not with the correct coefficient. 

Although Eq. (99) looks similar to the ghost and gluon loop contributions 
(43) and (44), there is a crucial difference between the genuine one-loop inte- 
grals and the effective loop diagram. The latter contains an additional factor 
Nc.afiFji{...)G\{...) which behaves as a running coupling a{y) in the ultraviolet 
and ensures that the loop is logarithmically subleading. 

The effective stmset contribution is constructed in exactly the same way as 
depicted above for the squint diagram. Note that the coefficients Cijm of the 
effective kernels for the sunset diagram are an order of magnitude smaller than 
those of the squint diagram, which agrees with the earlier conclusion that its 
infrared contribution is negligible. 

The introduction of the effective loops allows for a numerical study of the 
coupled ghost-gluon system including the main features of the two-loop diagrams 
to the gluon vacuum polarization. 

5.2 Results 

The existence of power laws satisfying the infrared consistency conditions does 
not yet guarantee that the coupled integral equations have solutions over the 
whole momentum region connecting the predicted infrared asymptotic behavior 
with the ultraviolet perturbative regime. A straightforward way to answer this 
question is to look numerically for solutions of the coupled equations. 
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To solve Eqs. (19) and (20) numerically the unknown seeds Z3 and of 
the equations are eliminated by subtracting each equation at a fixed subtraction 
point: 

' ■- n(p^) + n(a^) , 

(104) 



^ ^ E(p2) + E(4), 



and the renormalization conditions are fixed^ by requiring a{iJ,'^) = and 

GRiti-") = FR{^J?) = 1. 

The numerical method employed to solve the set of coupled integral equations 
(104) was developed in Ref. [37] and is described in more detail in Appendix F. 
Numerical solutions were sought for the ten truncations of the ghost-gluon system 
described in Table 3 of Sect. 4.4. As shown in that table, the lowest lying k of 
each truncation corresponds to a consistent numerical solution, while the higher 
value is spurious as no numerical solution having such an infrared behavior is 
found. The numerical solutions for the various parametrizations are very similar 
to one another, and typical SU(3) solutions^^ for truncation B.4 of Table 3 are 
shown in Figs. 8 and 9. The infrared power behaviors of the dressing functions 
and the infrared fixed point of the running coupling are in perfect agreement 
with the predictions of the infrared analysis, with k = 0.36 and ao = 2.4 for 
SU(3). Note that the infrared fixed point of the running coupling in Fig. 9 is the 
result of a product of power laws with vanishing total exponent (sec Eq. (52)), 
and its truthful rendering requires a numerical method with high accuracy. 

Because the equations are scale invariant, the momentum units of the so- 
lutions are arbitrary, and the scale should be fixed using experimental input. 
A typical choice is to fix the coupling at a specific perturbative mass scale, for 
instance a(M^). However, as this is a study of pure gauge theory and no quark 
loops are taken into account, the coupling does not run like the physical one, and 
such a scale matching would have no phenomenological use. Therefore, the mo- 
mentum units in Figs. 8 and 9 are left arbitrary for the time being. Nevertheless, 
it is clear that the solutions contain an intrinsic scale, which signals the tran- 
sition from the perturbative to the nonperturbative region and is fixed by the 
renormalization conditions. Note that because the equations satisfy multiplica- 
tive rcnormalizability the solutions for arbitrary renormalization conditions can 
easily be derived from the solutions computed for a particular renormalization. 

Of major interest is the comparison between recent SU(2) lattice simulations 
[17, 18, 19] and the Dyson- Schwinger results. The lattice data also exhibit in- 
frared power behaviors for the gluon and ghost dressing functions, and they seem 
to indicate that k = 0.5 as the gluon propagator goes to a non-zero, finite value 
at zero momentum. Furthermore, the lattice results feature a pronounced bump 



For reasons of numerical stability, the ghost equation is subtracted at ctq = 0, where 1/Gr{0) = 
0, and the renormalization conditions are chosen by fixing a,j, Fi{{a^) = Fa, and the coefHcicnt 
A of the asymptotic infrared power law for F{x). Because of multiplicative renormalizability 
these conditions are equivalent to the conditions a{n^) = afi,GR{iJ^) = Fr{ij,^) = 1. 
^"The relation between SU(N) and SU(3) results is trivial as can be seen from the equations of 
Sect. 2 and the discussion after Eq. (53). The solutions for F and G remain unchanged, while 
the running coupling a is multiplied by a factor 3/N. 
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X X 

Figure 8. Numerical solutions for truncation B.4 of Table 3 with k = 0.36. Left pane: Gluon 
dressing functions as function of squared momentum, with renormalization condition -F'(l) = 1. 
Right pane: Ghost dressing functions as function of squared momentum, with G(l) = 1. 
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Figure 9. Running coupling for trun- 
cation B.4 of Table 3, computed us- 
ing Eq. (35) and plotted as function 

of squared momentum. The infrared 
0.0001 0.01 1 100 10000 ^ , . , . 

nxed pomt value is ao = 2.4. 

X 

in the gluon dressing function, in the transition from the perturbative to the 
nonperturbative behavior. 

Because of its simphcity, the MR truncation cannot be expected to predict 
the correct value of the infrared exponent k. Moreover, for a given k the effective 
one-loop approximation to the two- loop diagrams (see Sect. 5.1) most probably 
affects the intermediate behavior of the solutions. An important goal of this in- 
vestigation is to find out to what extent the two-loop improved MR truncation of 
the ghost-gluon system is able to accommodate the various features exhibited by 
the lattice data. For this purpose the physically motivated parametrizations of 
the squint diagram proposed in Table 2 are now relaxed, and a, k = 0.5 solution 
is constructed by tuning the parameters. Note that the numerical solutions of 
the Dyson- Schwinger equations are stable with respect to variations of the pa- 
rameters. The main effects are changes in the infrared exponent and the related 
value of the infrared fixed point of the coupling, and changes in the steepness of 
the transition from perturbative to nonperturbative region, which is related to 
the height of the bump in the gluon dressing function. 

The MR truncation requires the specification of the ghost self-energy trun- 
cation as discussed in Sect. 4.1. As illustrated in Fig. 1 the value of the infrared 
fixed point of the coupling is related to the value of k through the ghost equa- 
tion conditions (60)- (62). For k = 0.5 the a-truncation (62) yields a higher value 
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ao = 5.2 for the SU(2) infrared fixed point (compared to 3.9 for the bare trun- 
cation), and seems in better agreement with the lattice data [19]. 

The squint diagram truncation is parametrized by the sets a4g and h4g of Eq. 
(46), and its infrared contribution is determined by k = K{2h4g—a4g). In practice, 
a K = 0.5 solution is constructed by solving the consistency equation (91) for k, 
after imposing the conditions fei + ^2 + ^3 + ^4 = —2k and k2 = ks (symmetry 
of the inner loop). The solution is not unique and one possible solution, which 
will be investigated in more detail, is k = (—2.072,0.536,0.536,0). 

The effective one-loop integrals (99) approximating the two-loop integrals in 
the numerical calculation, depend on k2 and ^3, which fix the coefficients Cijm 
of the effective integration kernel (101), and on ay, az, by and bz defined in Eq. 
(100). Given a set k, determining the infrared behavior of the squint diagram, 
there is one more parameter from the set Uy, az, by, that can be chosen 
freely. The remaining parameters are then fixed by Eqs. (102) and (103). This 
freedom of choice leaves the leading infrared and ultraviolet behaviors unchanged, 
but can influence the behavior of the solutions in the intermediate momentum 
region. Hence, it can be used as phenomenological parameter to recover some of 
the information that is unavoidably lost by replacing the two-loop integrals by 
effective integrals. In practice, the parameter by is tuned to reproduce the height 
of the bump in the lattice gluon dressing function, and good agreement is found 
for by = —0.6, which yields ay = 0.8, = 3.2 and bz = 1.6. 

The DSE results are compared with the SU(2) lattice data of ref. [19]. These 
lattice simulations were performed on a 16^ x 32 grid using 6 different /3-values 
ranging from 2.1 to 2.5. The momentum scale of the lattice calculations is fixed 
using a string tension of 440 MeV. The momentum units of the DSE results 
are fixed by maximizing the overlap between the DSE solution for the running 
coupling and the corresponding lattice data. The running coupling on the lattice 
was also computed from the dressing functions using Eq. (35) (see ref. [19]). 
Figures 10 and 11 compare the lattice data with the DSE results^^ computed 
with the two-loop improved MR truncation, and with the truncation of ref. [24] . 

The ability for the DSE formalism to reproduce the lattice data, by introduc- 
ing the two-loop improved MR truncation, is way beyond the initial expectations. 
Even though the DSE truncation was fixed using some characteristics of the lat- 
tice results, the DSE solutions are not molded to those on the lattice, and they 
genuinely predict a number of additional features which are consistent with the 
lattice data, as is discussed below. 

Although the height of the bump in the gluon dressing function (upper pane 
of Fig. 10) was tuned to that of the lattice, the position and width of the peak 
are genuine DSE predictions which agree well with the lattice data. The ghost 
dressing function (lower pane of Fig. 10) is in total agreement with the lattice 
data, and reproduces well the smooth evolution from the perturbative region 
into the 1/p infrared singularity. 

Even more remarkable is the agreement of the running coupling (Fig. 11) 
over the complete momentum region. The value of the fixed point is a DSE 



The data are now plotted on linear scale for both axes as the lattice data cover a much smaller 
momentum region than the DSE computations. Also note that the abscissa is p instead of p^. 
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Figure 10. Comparison of the lattice data from rcf. [19] with the DSE results using the two-loop 
improved MR truncation (with k = 0.5) and from rcf. [24] . Upper pane: Gluon dressing functions 
as function of momentum. Lower pane: Ghost dressing functions as function of momentum. 

prediction, determined by the ghost infrared condition (58). Even though the 
lattice simulations are not yet sampling the infrared region deep enough to clearly 
show the bend-over towards a fixed point, the existence and value of the DSE 
fixed point are consistent with the actual infrared behavior of the lattice data. In 
the DSE study, the steepness of a{x) in the transition region appears to be tightly 
related to the height of the bump in the gluon dressing. The computations show 
that the steepness of the DSE coupling agrees well with the lattice data precisely 
when the height of the gluon bump does too, and hence, their interrelation is 
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Figure 11. Running coupling as function of momentum for SU(2). Comparison of the lattice 
data from ref. [19] with the DSE results using the two-loop improved MR truncation (with 
K = 0.5) and from ref. [24]. The fixed point of the MR truncation is ao = 5.24. 

correctly reproduced by the DSE truncation. Tlie perturbative tail agrees very 
well with the lattice results too, illustrating that the solutions of the two-loop 
improved DSE truncation go beyond the leading logarithm. This contrasts with 
the usual one-loop inspired truncations, which typically yield ultraviolet tails 
that are 20% larger than the lattice results, as can be seen in Fig. 11 from the 
DSE curve taken from ref. [24]. 

It is also useful to construct a fit to the running coupling solution, which 
could be used in hadronic phenomenology: 

a{x) = — ^-^ [(co + cix'^) ao + x'^ a2e{x + r)] , (105) 

with X = P^/Aqqj-,, and a2i{x) is the two-loop perturbative running coupling 
[38]. The fixed point ao = 5.24 is known from the infrared analysis, r is set 
equal to 2, and the other parameters are fitted to the DSE solution of Fig. 11: 
Aqcd = 856 MeV, cq = 1.16, p = 0.66 and ci = -0.070. 

Note that other truncation parametrizations, which also give k, = 0.5 and re- 
produce the gluon bump, yield solutions which only differ in the nonperturbative 
approach to the fixed point of the running coupling, as shown in Fig. 12. 

The results presented in Figs. 10 and 11 are solutions for pure SU(2) gauge 
theory. For nature's QCD the corresponding SU(3) value of the fixed point is 
ao = 3.49. Furthermore, to get the correct running of the coupling, quark loops 
have to be included in the gluon vacuum polarization taking into account their 
dynamically generated mass. This is done by solving the coupled quark-ghost- 
gluon system in an MR truncation merging the truncation presented herein with 
those for the quark loop [25] and for the quark equation [34]. The solution of 
this coupled system will be reported elsewhere [39] , but from a previous study of 
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Figure 12. Running coupling for three 
different parametrizations, constructed 

to yield n = 0.5 and reproduce the 
bump in the gluon dressing function, 
and specified in the table below. 



a (-2.072,0.536,0.536,0) -0.6 
b (-1.552,0.526,0.526,-0.5) -0.15 
c (0.596,-0.5,-0.5,-0.596) -0.1 



the quark equation [34] it is already clear that a running coupling similar to Eq. 
(105) is capable of producing a strong enough breaking of the chiral symmetry 
to allow for consistent hadronic phenomenology. 



6 Conclusion 

In ref. [25] a truncation of the Landau gauge Dyson-Schwinger equations for 
the ghost and gluon propagators was proposed, which respects the principles of 
multiplicative renormalizability and reproduces the leading order resummed per- 
turbative results. This was achieved by directly approximating the integration 
kernels, rather than introducing explicit approximations to the full vertices and 
their corresponding renormalization constants. A remarkable feature of the trun- 
cation is that all diagrams in the gluon vacuum polarization contribute to the 
leading power in the infrared analysis. However, in that study the contributions 
of the two- loop gluon vacuum polarization diagrams could not be evaluated, and 
the MR truncation did not yield solutions for the coupled ghost-gluon system 
when only taking into account the one-loop diagrams. 

The present paper generalized further the MR truncation scheme of ref. [25] 
and confirmed that the interplay between all the loops in the gluon vacuum 
polarization is a key ingredient to find consistent MR solutions of the coupled 
ghost-gluon system. 

An analytic method was developed to compute the contributions of the two- 
loop diagrams to the leading infrared power behavior of the gluon vacuum polar- 
ization, such that a complete infrared analysis of the coupled ghost-gluon system 
could be performed. The investigation showed that the squint diagram is able to 
supply an infrared contribution that allows for propagator solutions with infrared 
power behavior. The gluon and ghost dressing functions respectively behave like 
(jP'Y'K g^j^j^ ^^2-j-K^ g^j^j ^YiQ running coupling exhibits an infrared fixed point. 
The value of the infrared exponent depends on the details of the truncation, and 
some simple realizations of the truncation yield values of k G [0.17, 0.53]. 

The numerical study of the coupled ghost-gluon system was facilitated by 
approximating the two-loop gluon vacuum polarization by an effective one-loop 
integral with the same leading infrared behavior. The numerical analysis pro- 
duced solutions of the ghost-gluon system over the whole momentum range, and 
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their infrared behaviors agree with the predictions of the analytic infrared anal- 
ysis. 

Recent lattice simulations suggest that k = 0.5, and DSE solutions having 
such an infrared behavior were constructed by adjusting the truncation param- 
eters. The numerical solutions for the dressing functions and running coupling 
showed that the two-loop improved MR truncation of the Dyson-Schwinger equa- 
tions is able to produce solutions that are in extremely good agreement with the 
lattice data. 

This paper presented a study of the pure gauge sector. In a forthcoming 
study [39], this MR truncation will be merged with that of the quark equation 
developed in ref . [34] , to allow for a complete self-consistent study of dynamical 
chiral symmetry breaking in QCD. 
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Appendix A: Color Contraction 

The color contraction identities needed in Sect. 2 are 

fciC2afciC2b — CA^ab , (-^•-^) 
/ciC2a/ciC36/c2C3C = ~^fabc , (■^•2) 

cl 

faciC2 fbc2C4 fciC3C5 fc2C4,C5 — ~^^ab • (-^•^) 

Appendix B: Multiplicative Renormalizability in Quenched QED 

To get some feeling for the truncation suggested in Sect. 3, it is good to make a detour over 
quenched QED [29] and see how there the multiplicative renormalizability of the massless quark 
propagator is achieved in the Dyson-Schwinger equation, at least on the perturbative level. 

Although the divergence in the vertex can be written as a multiplication of the bare vertex 
by a divergent product of dressing functions, a variety of studies have shown that in general 
such a dressed vertex is not sufficient in its application to the Dyson-Schwinger equations to 
retrieve the correct divergence and the multiplicative renormalizability of the propagators. 

The DSE for the quark dressing function is: 

Expanding the RHS to order a, i.e. replacing the full vertex by its bare analog, yields 

' -i-^m^. (B.2) 



Zip'^) 47r A2 

To 0{a) this gives 

47r A2 



ZiP)-^ + ^^n^, (B.3) 



and multiplicative renormalizability tells us that 



Z{P) = ( ) , (B.4) 
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after resummation of the leading- logarithm contributions. 

The question now is: how is the correct In^-term produced such that the series (B.3) is 
continued in a way consistent with the MR solution (B.4)? This was discussed previously in 
ref. [29], and the mechanism is clarified in detail below. 

For large momentum the leading- logarithm one- loop fermion-photon vertex is [29] : 



7m 



1 - 



7^ In — 



a 
An 



1 

fc2 



In 



(B.5) 



Note that the ultraviolet divergent part of the vertex (A-dependent) is exclusively contained 
in the term proportional to 'y^ of the vertex, and for large k the vertex correction is given by 
1/Z{k^). Nevertheless, the divergent part of the vertex is not sufficient to generate the correct 
divergence of the propagator in the DSE, nor is it capable to produce an MR solution: the 
finite, momentum dependent part of the vertex plays a crucial role. 

The next-to-leading order term of Z(jp) is generated by substituting the 0{a) dressing 
function (B.3) and vertex (B.5) (neglecting subleading parts independent of ^) in Eq. (B.l). 
After performing the angular integrals, the ultraviolet integral becomes 



Z(p2) ^ 4n 



I" 



dk^ 
fc2 



+ 



dk' 

fc2 



-In 



A2 



+ ln 



k' 



+ In 



hi 

A2 



+ • 



(B.6) 



r(i) z(o) 



r(o)z(i) 



where r''^-', Z*"', r'^', Z'^' respectively correspond to the bare and one-loop corrections to the 
vertex and dressing function. Integration yields the dressing function to two-loop: 



Z(p2) 



= 1 - 



^In^ 
An A2 



4n A2 



(B.7) 



These are the first three terms of the Taylor expansion of the inverse of Eq. (B.4). From Eq. 
(B.6) it is clear that only taking into account the divergent parts of the one-loop vertex and 
dressing function is not sufficient to reproduce the correct two-loop contribution to the dressing 
function, which is crucial to ensure the resummability of the series required by multiplicative 
renormalizability. Therefore, DSE truncations where the full vertex is simply approximated 
by its divergent part, i.e. bare vertex multiplied with dressing functions, will not preserve the 
multiplicative renormalizability of the propagators. 

How can MR be preserved and the correct divergence be reproduced in a truncation using 
only the part of the vertex? 

The self-energy integral of Eq. (B.l) can be written as 



K{x,y,z)Z{y)= / [Z'\y)K{x,y, z)] 



1 



'Z{y) 



(a) 

where Z^K and 1/Z were expanded to 0{a). The interesting observation is that the factor 
(a) yields a vanishing In^ integral contribution. This explicitly shows that the In^-correction 
to 1/Z{x) is correctly generated by the integral over the bare vertex kernel Ko times the In- 
correction of 1/Z{y), i.e. 



/. 



K{x,y,z)Z{y) 



-L 



Ko{x,y.z)^^ 



to 0{\n) . 



(B.9) 



In fact, using the Taylor expansion of the power law (B.4) it is easy to see that the RHS of 
Eq. (B.9) reproduces term by term the leading-logarithm scries for 1/Z{x), and Eq. (B.9) is valid 
at each order in perturbation theory. This means that the MR truncation correctly reproduces 
the resummed leading- logarithm results. Note that this does not imply that the kernels Z^K 
and Ko are identical, but just that they yield the same ultraviolet integral contributions. In 
the same way, the MR truncation can be shown to reproduce the resummed leading-logarithm 
results of the gluon and ghost propagators in QCD. 
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The MR truncation of Sect. 3 assumes that the use of Eq. (B.9) can reasonably be extended 
to all external momenta, at least in a slightly more general form, as it ensures that multiplicative 
renormalizability is satisfied over the complete momentum domain. 

Appendix C: Feynman Integration Trick 

To solve the infrared integrals of Sect. 4 using the integrals given in Appendix D, it is helpful to 
transform denominator products of squared momenta in the kernels using the following equality: 

[\te-\i-tr' ^ (c.i) 

qr^qr' . . . - T{m,)T{m,) . . . r(m^ j, ■ ■ ■ dtr. 5 - l) ^^^^^^ . 

(C.2) 

Eventually the Feynman parameter t is integrated out using the well-known formula: 

\tt%l-tf=^±±^^. (C.3) 



p2-q2b r{a)r{b)Jo ^ ' [tp2 + (1 _ f)g2] 

which can be generalized as 



Appendix D: Integrals 

The tensor integrals needed to compute the infrared contributions of the one- and two-loop 
diagrams are [40] 

Jin, A) ^ I j,^^ = ^iA^^A-", (D.l) 
A) = / = -^-^A-"% , (D.2) 



1 


1 


r(n - 2) 


{P + A)" 


167r2 


r(n) ■ 




1 


r(n - 3) 


(£2 -1- A)" 


167r2 


r(n) ■ 


p f e e 


1 


r(n - 4) 




167r2 


r(n) • 






i r(n - 



J^.pa(n, A) = ^ {p+AY " 16^2 ^^^FH^^" " J(<5m-V + ^i^pSu. + 5^^5^p) , (D.3) 

~i~ pu^^ij ^^ ^p^f^vcr ~\~ ^ pa ^vy^^T p i^^pu^p^f ^ ^pp^u^ ~\~ ^ p'-^^u p)^T(T^ , 

and the odd-rank tensor integrals vanish because of symmetry properties. 

By applying the Feynman integration trick (C.I) and using the integration formula (D.l), 
a typical scalar integral, often needed in one-loop integrals, is given by 

T(n M - / ^ -±- - 1 r(2-a)r(2-b)r(a + b-2) 

^ ' ^ ~ y (277)4 2/";?'' - 167r2 r(a)r(6)r(4 - a - 6) ' ^ ' 

with X = ., y = q^ and z = {p — q)^. 

Appendix E: Effective Two-loop Kernels 

As discussed in Sect. 5.1 the two-loop diagrams are approximated by effective one-loop integrals 
in the numerical calculation. The kernel of the effective squint diagram can be written as 



12_3y_4 _C!f; + liLl -I- -I- 20^ 

Qfix,y,z)=A{k2,k3){ ^ ^ + 36 + ^+ ^ + 



(E.l) 

with 

-r(l + fca) r(l fca) r(-fc2 - fca) (-5 3 fca' 3 fca' fca fcg fca fca) 

■^^'^ 256 7r2r(2-fe)r(2-ife3)r(3 + fc2 + fc3) ■ ^ ' 
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Kernel (E.l) contains potential quartic (l/x^-term) and quadratic (l/x-term) divergences. 
Although these terms formally cancel after integration, their numerical implementation can 
yield serious problems because incomplete cancellations will give rise to residual unrenormal- 
izable divergences [37, 41]. A simple, elegant solution consists in writing the expression over a 
common denominator and then express z as x + y — 2y/xycoa0: 



Qf(x,y,z)=A(k2,k3) ^ (-1 + 4cos^ 6») - 20y^/iy cos6> (l + 8cos^ 6i) 

xyz^ ^ ^ ^ ' (E.3) 

+ 4a;j/(-7 + 43cos^6») -36 Xyfxy cos Q ] . 

Once written in this form the simple trigonometric form of the remaining l/x-term ensures that 
the quadratic divergence will be cancelled properly. An analogous expression can be written for 
the sunset diagram. 

Appendix F: Numerical Method 

The numerical method employed to solve the set of coupled integral equations (104) was first 
devised for the study of dynamical chiral symmetry breaking in strong couphng QED [37] . After 
changing to spherical coordinates and integrating out two angles, the integral equations can be 
written as (dropping the subscript R on the renormalized dressing functions) 



1- = ^+ Tdy r deKirlG,F]{x,y,e), 

^[X) Jo Jo ,p - 

-1 1 POO pTT \ ■ / 

7T^ = 7^+ dy de Ka[G,F]{x,y,e). 

G(x) Gct Jo Jo 



The unknown dressing functions are approximated using Chebyshev expansions [42] in 
log(a;) for x € [e^, A^]: 

N-l 

logF(x) :=|^+^a, T,(.(x)), 

(F.2) 

logG(x) :=| + ^&,T,(s(x)), 

where and A^ are infrared and ultraviolet parameters introduced for numerical purposes only, 
and 

log(x-/A6) 

log(A/e) ' ^^-^^ 
such that s £ [—1,-1-1] when x £ [e^,A^]. The number of Chebyshev polynomials is chosen 
between 30 and 60, depending on the steepness of the running coupling in the intermediate 
momentum regime and the required accuracy. In the infrared, for x < , the dressing functions 
are taken to obey their asymptotic power law behaviors 

F{x) ■.= Ax^\ 

(F.4) 

G{x):= Bx-" , 

and in the ultraviolet, for a; > A^, the functions are frozen to their values at A^ given by the 

Chebyshev expansions (F.2). 

The system of equations (F.l) is required to be satisfied at A'^ external momenta 



= Ae I — ) , Si = cos 



(i - l/2)7r 



N 



i = l,...,N, (F.5) 



where the Si are chosen to be the zeroes of the Chebyshev polynomial Tjv of degree N, and the 
integrals are approximated by suitable numerical quadratures: 

■pf-T = ■pr + Y^^WjWkKF[F,G\{xi,yj,ek), 
q7^ = ^ + XIX] Wjw;j^'f^G[F,G] (a;i,j/j, 61/0 ), 

^ ^ j=l k=l 
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with radial integration points yj, j = l,...,Ny, and angular integration points 9k, k = 
1, . . . , Ne- Note that the nodes yj can depend on the external momentum Xi. 

The radial integrals are split into three parts: [0, e^], [e^,A^], and [A^^,oo]. An accurate 
evaluation of the infrared integrals [0, e^] is important to reproduce the infrared power laws, 
and they are solved using a Gauss-Legendre integration with 64 points. The major integration 
range [e^, A^] has to be handled caxefuUy as the integrand is not smooth at y = x. A composite 
Gauss-Legendre rule is used (after transformation to t = \og{y)), with main nodes located at 
e^, and the A'^ values of external momenta, and an M-point Gauss-Legendre rule is used 
between adjacent nodes. The integrand is concentrated around y = x, and hence the number of 
points in the elementary rule is adapted accordingly: a high accuracy rule with M = 32 is used 
when the integration frame overlaps with the interval [x/10, lOx], while otherwise M = 4. For 
the subtracted equations the ultraviolet integral [A^, oo] is negligible. The angular integrals are 
computed using a A^g-point Gauss-Logondro integration, with Ng — 33. 

The method described above transforms the set of coupled integral equations (F.l), into 
a set of 2N non-linear algebraic equations (F.6) for the Chebyshev coefficients aj,bj of the 
dressing functions. Symbolically this can be written as 



f{x)[F,G]=Q 
g{x)[F,G] = 



fi[aj,bj 







gi[aj,bj] = 



(F.7) 



The latter set of equations is solved using the Newton iterative method, which is quadrat- 
ically convergent for starting guesses close enough to the solutions. Successive iterations are 
given by 



"3 - 

6"+' = 6" 



A6"+ 



(F.8) 



where the increments Aaj and Abj are solutions of the set of 2N linear equations 



dgi [a 



An- 



■Ab] 



-1 _|_ dgi[a 



dbi 



Abj^ =gi[aj,bj]. 



i = l. 



.AT, 



(F.9) 



and summation over the repeated index j is assumed. The quadratic convergence property is 
clearly observed in the numerical implementation as the distance between successive iterations 
decreeises by an order of magnitude at each iteration step in the vicinity of the solution. Hence, 
the method is fast and gives a reliable estimate of the accuracy of the solutions. 

Note that the quadratic convergence only sets in "close" enough to the exact solution, and 
good initial guesses dramatically speed up the convergence rate. These can easily be constructed, 
as the infrared and ultraviolet asymptotic behaviors of the dressing functions are known. Typical 
starting guesses are 

2k 



Fix) 



G{x) = B 



x/xt + 1 



x/xt + 1 



with 



a{t) 



co + t^ 



Co ao + 



47r 



&(x/xt) 
a{x/xt 



1 

logt 



13/22 



9/44 



(F.IO) 



1 - t 



where co and Xt can be tuned to improve the starting guess. The procedure (F.8) is iter- 
ated till the convergence criterion maxj (| Ao"+^|) < 10~^ maxj(|a"+^|) and maxj(|Afe"+^|) < 
10~^ maxj(|6^"'"'^|) is satisfied. Starting from pure power laws (F.4) for all momenta typically 
requires about 20 iterations to converge, while the more refined starting guesses (F.IO) reduce 
the iteration procedure to 4 iteration steps. 

The total accuracy of the method is given by a combination of the accuracies of the Cheby- 
shev expansion, the quadrature rules, and the Newton iteration method. 
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